'''
实现影像重投影和拼接
'''
import datetime
import glob
import os

from osgeo import gdal, gdalconst
from tqdm import tqdm


# 实现合并函数
def merge_clip(merge_dir1, merge_dir2, output_dir):
    # 获取待合并文件夹下的所有tif文件
    merge_file1 = merge_dir1 + "/*.tif"
    merge_file2 = merge_dir2 + "/*.tif"
    # 获取文件夹下所有tif文件
    file_list1 = glob.glob(merge_file1)
    file_list2 = glob.glob(merge_file2)
    dataset = gdal.Open(file_list1[0])
    project = dataset.GetProjection()
    nodata = dataset.GetRasterBand(1).GetNoDataValue()
    # print(file_list1)
    for i in tqdm(range(len(file_list1))):
        output_tif = output_dir + "/" + file_list1[i].split("\\")[-1].split("_")[0].split('.')[1][1:] + ".tif"
        # gdal.BuildVRT(output_tif, [file_list1[i], file_list2[i]])
        gdal.Warp(output_tif, [file_list1[i], file_list2[i]],
                  srcSRS=project,
                  dstSRS="EPSG:4326",
                  format="GTiff",
                  # outputBounds=[91, 31, 94, 34],
                  outputBounds=[98, 38, 101, 42],
                  srcNodata=nodata,
                  dstNodata=nodata,
                  # cutlineDSName=r"G:\MOD13A2\shp\study_area.shp",
                  # cropToCutline=True,
                  resampleAlg=gdalconst.GRA_NearestNeighbour)
    print("merge-clip success")


# 实现合并函数
def merge(merge_dir1, merge_dir2, output_dir):
    # 获取待合并文件夹下的所有tif文件
    merge_file1 = merge_dir1 + "/*.tif"
    merge_file2 = merge_dir2 + "/*.tif"
    # 获取文件夹下所有tif文件
    file_list1 = glob.glob(merge_file1)
    file_list2 = glob.glob(merge_file2)

    dataset = gdal.Open(file_list1[0])
    project = dataset.GetProjection()
    nodata = dataset.GetRasterBand(1).GetNoDataValue()
    print(nodata)
    for i in tqdm(range(len(file_list1))):
        # dataset1 = gdal.Open(file_list1[0])
        # dataset2 = gdal.Open(file_list2[0])
        output_tif = output_dir + "/" + file_list1[i].split("\\")[-1].split("_")[0].split('.')[1][1:] + ".tif"
        gdal.Warp(output_tif, [file_list1[i], file_list2[i]],
                  srcSRS=project,
                  dstSRS='EPSG:4326',
                  format="GTiff",
                  srcNodata=nodata,
                  dstNodata=nodata,
                  resampleAlg=gdalconst.GRA_NearestNeighbour)
    print("merge success")


def clip(input_dir, output_dir):
    # 获取待合并文件夹下的所有tif文件
    input_file = input_dir + "/*.tif"
    # 获取文件夹下所有tif文件
    file_list = glob.glob(input_file)
    dataset = gdal.Open(file_list[0])
    nodata = dataset.GetRasterBand(1).GetNoDataValue()
    # print(file_list1)
    for i in tqdm(range(len(file_list))):
        # output_tif = output_dir + "/" + file_list[i].split("\\")[-1].split('.')[0] + ".tif"
        # output_tif = output_dir + "/" + file_list[i].split("\\")[-1].split("_")[0].split('.')[1][1:] + ".tif"
        date = os.path.basename(file_list[i]).split(".")[0]
        # 根据日期计算儒略日
        # targetday = datetime.date(int(str(date)[0:4]), int(str(date)[4:6]), int(str(date)[6:8]))
        # day = targetday - datetime.date(targetday.year - 1, 12, 31)
        # julian_day = str(date)[0:4] + str(day.days)
        output_tif = output_dir + "/" + date + ".tif"
        gdal.Warp(output_tif, file_list[i],
                  dstSRS="EPSG:4326",
                  format="GTiff",
                  # outputBounds=[91, 31, 94, 34],
                  # targetAlignedPixels=True,
                  outputBounds=[145.5, -35.5, 148, -34.5],
                  # xRes=0.25,
                  # yRes=0.25,
                  srcNodata=nodata,
                  dstNodata=nodata,
                  resampleAlg=gdalconst.GRA_NearestNeighbour)
    print("clip success")


if __name__ == '__main__':
    # LST合并和裁剪
    # merge_dir1 = r"G:\test\process_result\LST_tif_new\h25v04\LST_day"
    # merge_dir2 = r"G:\test\process_result\LST_tif_new\h25v05\LST_day"
    # output_dir = r"G:\test\process_result\LST_tif_new\merge_clip\LST_day"
    # if not os.path.exists(output_dir):
    #     os.makedirs(output_dir)
    # merge_clip(merge_dir1, merge_dir2, output_dir)
    # merge(merge_dir1, merge_dir2, output_dir)
    # clip_input_dir = r"G:\test\process_result\LST_tif_new\merge\LST_day"
    # clip_output_dir = r"G:\test\process_result\LST_tif_new\clip\LST_day"
    # if not os.path.exists(clip_output_dir):
    #     os.makedirs(clip_output_dir)
    # clip(clip_input_dir, clip_output_dir)

    # 裁剪
    # clip_dir = r"G:\test\process_result\LST_rec_tif\aqua\merge\night"
    # output_dir = r"G:\test\process_result\LST_rec_tif\aqua\clip\night"
    # clip(clip_dir, output_dir)

    # EVI合并和裁剪
    # merge_dir1 = r"G:\test\process_result\EVI_tif\h25v05\QC"
    # merge_dir2 = r"G:\test\process_result\EVI_tif\h26v05\QC"
    # output_dir = r"G:\test\process_result\EVI_tif\merge_clip\QC"
    # if not os.path.exists(output_dir):
    #     os.makedirs(output_dir)
    # merge_clip(merge_dir1, merge_dir2, output_dir)

    # GLASS合并和裁剪
    # merge_dir1 = r"G:\test\process_result\GLASS_tif\h25v05\QC_VIS"
    # merge_dir2 = r"G:\test\process_result\GLASS_tif\h26v05\QC_VIS"
    # output_dir = r"G:\test\process_result\GLASS_tif\merge_clip\QC_VIS"
    # if not os.path.exists(output_dir):
    #     os.makedirs(output_dir)
    # merge_clip(merge_dir1, merge_dir2, output_dir)

    # CCI裁剪
    # clip_dir = r"G:\test\SMAP_NEW\ESACCI\TIF\2015-2017"
    # output_dir = r"G:\test\SMAP_NEW\ESACCI\TIF\2015-2017-CLIP"
    # if not os.path.exists(output_dir):
    #     os.makedirs(output_dir)
    # clip(clip_dir, output_dir)

    input_path = r"G:\test\SMAP_NEW\svtm_nsw_extant_pct_vc1_1_m1_1\SVTM_NSW_Extant_vC1_1_M1_1_5m.tif"
    output_path = r"G:\test\SMAP_NEW\svtm_nsw_extant_pct_vc1_1_m1_1\SVTM_NSW_Extant_vC1_1_M1_1_5m_clip.tif"
    input_dataset = gdal.Open(input_path)
    geo = input_dataset.GetGeoTransform()
    proj = input_dataset.GetProjection()
    nodata = input_dataset.GetRasterBand(1).GetNoDataValue()
    gdal.Warp(output_path, input_path,
              srcSRS=proj,
              dstSRS='EPSG:4326',
              outputBounds=[145.5, -35.5, 148, -34.5],
              format="GTiff",
              srcNodata=nodata,
              dstNodata=nodata,
              resampleAlg=gdalconst.GRA_NearestNeighbour)


